Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_CaboVerde.50	Actinomycetota	Actinomycetia	Actinomycetales	Actinomycetaceae	NA	0.953
bin_Aruba.25	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium pseudocatenulatum	0.995
bin_Aruba.2	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium adolescentis	0.994
bin_Aruba.6	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium longum	0.985
bin_Denmark.14	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium gallinarum	0.950
bin_Aruba.47	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	NA	0.964
bin_Aruba.57	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	Corynebacterium pyruviciproducens	0.931
bin_Aruba.51	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	Parolsenella sp900544655	0.969
bin_Denmark.44	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	UBA7748 sp900314535	0.951
bin_Aruba.14	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Denmark.4	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Aruba.21	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	NA	0.999
bin_Brazil.61	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella tanakaei	0.998
bin_Aruba.39	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp902362275	0.992
bin_Spain.84	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900555555	0.992
bin_Brazil.151	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella intestinalis	0.991
bin_Denmark.127	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900548365	0.991
bin_Denmark.17	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella bouchesdurhonensis	0.983
bin_Brazil.81	Bacillota	Bacilli	RFN20	CAG-826	NA	0.907
bin_Brazil.109	Bacillota_A	Clostridia	Oscillospirales	Butyricicoccaceae	AM07-15 sp003477405	0.936
bin_Aruba.28	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.995
bin_Malaysia.9	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Dysosmobacter welbionis	0.986
bin_Malaysia.116	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Flavonifractor plautii	0.981
bin_Malaysia.34	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Lawsonibacter sp000177015	0.954
bin_Aruba.54	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.959
bin_CaboVerde.35	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.930
bin_Aruba.36	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	Peptoniphilus_C sp902363535	0.913
bin_Aruba.52	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.906
bin_Malaysia.26	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	NA	0.976
bin_Aruba.29	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Gemmiger variabilis_C	0.941
bin_Malaysia.103	Bacillota_C	Negativicutes	Acidaminococcales	Acidaminococcaceae	Acidaminococcus intestini	0.942
bin_Malaysia.8	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Dialister sp900543165	0.981
bin_Denmark.60	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Allisonella histaminiformans	0.938
bin_Aruba.64	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.970
bin_Brazil.58	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera elsdenii	0.969
bin_CaboVerde.38	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.961
bin_Malaysia.81	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera stantonii	0.943
bin_Malaysia.96	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Caecibacter sp003467125	0.919
bin_Brazil.163	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella lascolaii	0.992
bin_CaboVerde.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella copri	0.981
bin_Aruba.10	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900544825	0.978
bin_Brazil.5	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900540415	0.929
bin_Spain.21	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprophilus	0.929
bin_Malaysia.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp000437675	0.924
bin_Malaysia.117	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp900541335	0.919
bin_Malaysia.77	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900542985	0.913
bin_Brazil.75	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.964
bin_Malaysia.61	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	Helicobacter_C labetoulli	0.964
bin_Brazil.128	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.963
bin_Brazil.70	Desulfobacterota	Desulfovibrionia	Desulfovibrionales	Desulfovibrionaceae	Mailhella sp003150275	0.944
bin_Brazil.146	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp001917125	0.969
bin_Brazil.9	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp000436375	0.919
bin_Aruba.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.963
bin_Brazil.51	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Sutterella wadsworthensis_A	0.963
bin_Spain.43	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.959
bin_Brazil.64	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.948
bin_Brazil.92	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.948
bin_Brazil.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.930
bin_CaboVerde.33	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum sp900543125	0.995
bin_Brazil.82	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum succiniciproducens	0.990
bin_CaboVerde.55	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum_A thomasii	0.939

8.3.2 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_Denmark.27	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_B hirae	0.000
bin_Brazil.12	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus faecalis	0.001
bin_Brazil.170	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	NA	0.014
bin_CaboVerde.24	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_E cecorum	0.020
bin_CaboVerde.16	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_Denmark.56	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Levilactobacillus brevis	0.000
bin_Malaysia.72	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_CaboVerde.60	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Malaysia.127	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Denmark.61	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Latilactobacillus sakei	0.003
bin_Malaysia.35	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus animalis	0.003
bin_Malaysia.4	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus ruminis	0.035
bin_CaboVerde.14	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus acidophilus	0.067
bin_Aruba.18	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus lutetiensis	0.001
bin_Brazil.22	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus pasteurianus	0.001
bin_Denmark.117	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Lactococcus garvieae	0.001
bin_Denmark.113	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus alactolyticus	0.002
bin_Brazil.69	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Clostridium_A leptum	0.097
bin_Denmark.93	Bacillota_A	Clostridia	Lachnospirales	Anaerotignaceae	Anaerotignum sp001304995	0.074
bin_CaboVerde.3	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Denmark.42	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Brazil.136	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium thermobutyricum	0.041
bin_Brazil.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900555395	0.026
bin_Brazil.89	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900550235	0.034
bin_Denmark.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerostipes hadrus	0.040
bin_Denmark.19	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900556835	0.052
bin_Denmark.118	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_A longicatena	0.067
bin_Malaysia.108	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.081
bin_Brazil.113	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerobutyricum sp900754855	0.095
bin_Denmark.50	Bacillota_A	Clostridia	Peptostreptococcales	Peptostreptococcaceae	Peptacetobacter sp900539645	0.053
bin_Brazil.107	Bacillota_A	Clostridia	Monoglobales_A	UBA1381	CAG-41 sp900066215	0.016
bin_Brazil.159	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900541465	0.040
bin_Brazil.140	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_A sp900543175	0.046
bin_Malaysia.6	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900545035	0.088

8.3.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

Rows: 202 Columns: 2
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): phylum, colors

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.5 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")

# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.5.0.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.5.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0103	0.008544525	0.0004896855	0.016616286	0.916	0.084	Biosynthesis	Nucleic acid biosynthesis	UDP/UTP
D0504	0.004425848	0.0003985114	0.009383968	0.938	0.062	Degradation	Amino acid degradation	Methionine
D0906	0.003599450	0.0001480281	0.007752579	0.926	0.074	Degradation	Antibiotic degradation	Fosfomycin
D0205	0.012498142	0.0019329087	0.023496427	0.944	0.056	Degradation	Polysaccharide degradation	Pectin
D0208	0.010205394	0.0017412487	0.019823177	0.938	0.062	Degradation	Polysaccharide degradation	Mixed-Linkage glucans

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.003892460	-0.009445156	-1.598674e-04	0.047	0.953	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021086408	-0.039777480	-3.107959e-03	0.077	0.923	Biosynthesis	Amino acid biosynthesis	Glutamate
B0302	-0.004702329	-0.010402985	-4.605576e-04	0.036	0.964	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.013024436	-0.024370029	-3.149799e-03	0.048	0.952	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011613334	-0.021373921	-2.413267e-03	0.064	0.936	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0804	-0.016684242	-0.031511256	-3.389357e-03	0.052	0.948	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016440795	-0.031679497	-2.316252e-03	0.070	0.930	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009340160	-0.018754997	-7.868225e-04	0.084	0.916	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.012342178	-0.023948940	-4.468990e-04	0.094	0.906	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002061541	-0.003660689	-5.929844e-04	0.040	0.960	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004458863	-0.007873548	-1.198957e-03	0.030	0.970	Degradation	Amino acid degradation	Ornithine
D0508	-0.003160979	-0.007357970	-1.369150e-04	0.088	0.912	Degradation	Amino acid degradation	Lysine
D0903	-0.003870200	-0.009394756	-1.560699e-04	0.047	0.953	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015263907	-0.027752320	-3.008654e-03	0.071	0.929	Degradation	Antibiotic degradation	Macrolide
D0601	-0.008936709	-0.017382607	-1.878995e-03	0.042	0.958	Degradation	Nitrogen compound degradation	Nitrate
D0611	-0.003870200	-0.009394756	-1.560699e-04	0.047	0.953	Degradation	Nitrogen compound degradation	Phenylethylamine
D0603	-0.001894066	-0.003700352	-3.047402e-04	0.053	0.947	Degradation	Nitrogen compound degradation	Urate
D0610	-0.002980443	-0.004927790	-8.438420e-04	0.059	0.941	Degradation	Nitrogen compound degradation	Methylamine
D0606	-0.005546897	-0.011162503	-6.350907e-04	0.073	0.927	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001551679	-0.002998413	-8.337605e-05	0.089	0.911	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001834932	-0.002282080	-9.835430e-05	0.009	0.991	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001834932	-0.002282080	-9.835430e-05	0.009	0.991	Degradation	Xenobiotic degradation	Xylene
D0807	-0.003989919	-0.008479592	-6.243554e-04	0.057	0.943	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004699653	-0.010697649	-5.137227e-04	0.060	0.940	Degradation	Xenobiotic degradation	Trans-cinnamate

8.5.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.003892460	-0.009445156	-1.598674e-04	0.047	0.953	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021086408	-0.039777480	-3.107959e-03	0.077	0.923	Biosynthesis	Amino acid biosynthesis	Glutamate
B0302	-0.004702329	-0.010402985	-4.605576e-04	0.036	0.964	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.013024436	-0.024370029	-3.149799e-03	0.048	0.952	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011613334	-0.021373921	-2.413267e-03	0.064	0.936	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0804	-0.016684242	-0.031511256	-3.389357e-03	0.052	0.948	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016440795	-0.031679497	-2.316252e-03	0.070	0.930	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009340160	-0.018754997	-7.868225e-04	0.084	0.916	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.012342178	-0.023948940	-4.468990e-04	0.094	0.906	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002061541	-0.003660689	-5.929844e-04	0.040	0.960	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004458863	-0.007873548	-1.198957e-03	0.030	0.970	Degradation	Amino acid degradation	Ornithine
D0508	-0.003160979	-0.007357970	-1.369150e-04	0.088	0.912	Degradation	Amino acid degradation	Lysine
D0903	-0.003870200	-0.009394756	-1.560699e-04	0.047	0.953	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015263907	-0.027752320	-3.008654e-03	0.071	0.929	Degradation	Antibiotic degradation	Macrolide
D0601	-0.008936709	-0.017382607	-1.878995e-03	0.042	0.958	Degradation	Nitrogen compound degradation	Nitrate
D0611	-0.003870200	-0.009394756	-1.560699e-04	0.047	0.953	Degradation	Nitrogen compound degradation	Phenylethylamine
D0603	-0.001894066	-0.003700352	-3.047402e-04	0.053	0.947	Degradation	Nitrogen compound degradation	Urate
D0610	-0.002980443	-0.004927790	-8.438420e-04	0.059	0.941	Degradation	Nitrogen compound degradation	Methylamine
D0606	-0.005546897	-0.011162503	-6.350907e-04	0.073	0.927	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001551679	-0.002998413	-8.337605e-05	0.089	0.911	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001834932	-0.002282080	-9.835430e-05	0.009	0.991	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001834932	-0.002282080	-9.835430e-05	0.009	0.991	Degradation	Xenobiotic degradation	Xylene
D0807	-0.003989919	-0.008479592	-6.243554e-04	0.057	0.943	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004699653	-0.010697649	-5.137227e-04	0.060	0.940	Degradation	Xenobiotic degradation	Trans-cinnamate

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.5.2.1 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D02	8.144008e-03	-0.003585743	0.0201994472	0.820	0.180
B08	7.729591e-03	-0.003619195	0.0188190890	0.782	0.218
B01	1.079193e-03	-0.006174439	0.0083360499	0.643	0.357
S01	4.929735e-04	-0.012595581	0.0138898138	0.556	0.444
D01	4.722996e-05	-0.005072307	0.0047107372	0.430	0.570
B09	8.550675e-05	-0.000537705	0.0005261702	0.375	0.625

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D08	-1.128130e-03	-0.0021383271	-0.0001797254	0.046	0.954
B03	-1.050395e-02	-0.0186128263	-0.0028600447	0.056	0.944
D06	-3.041585e-03	-0.0069243804	0.0001474005	0.112	0.888
B04	-8.359605e-03	-0.0175260241	0.0009514529	0.133	0.867
B06	-6.693422e-03	-0.0163287939	0.0032004529	0.184	0.816
D07	-1.184844e-02	-0.0305629696	0.0040330573	0.194	0.806
D05	-1.464921e-03	-0.0078777280	0.0038694072	0.214	0.786
D03	-3.992202e-03	-0.0126804243	0.0032883113	0.224	0.776
S03	-9.018483e-03	-0.0330709641	0.0152032371	0.250	0.750
B02	-3.531226e-03	-0.0122475416	0.0047519621	0.270	0.730
D09	-1.973050e-03	-0.0085424585	0.0046762375	0.307	0.693
S02	-4.605229e-03	-0.0136699132	0.0030797554	0.323	0.677
B07	-3.729566e-03	-0.0158138065	0.0091826160	0.330	0.670
B10	-8.894201e-06	-0.0003238718	0.0002411792	0.484	0.516

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")